/*
 * fieldTest.cpp
 *
 *  Created on: Jul 31, 2010
 *      Author: nish
 */
#include <iostream>
/*
 * problem - measure the potential profile of over 100 atoms laid in 10 by 10 matrix in a square,
 * shoot an electron with an initial kinetic energy just a little above the peak of the nuclear potential.
 * How is the monte carlo equation called semiclassical, as the force and velocity effects are treated as
 * newtonian.
 * Ultimate goal -> simulation of electron spin transport in quantum dots -> using quantum montecarlo
 * so -> we have electron spin transport model
 * quantum dot model -> mostly potential effects and collision/momentum transfer based effects.
 * quantum montecarlo -> transport model -> need to incorporate the electron spin transfer.
 * pit falls -> when reading material, can't understand much. i always need momentum, can't get stuck for long.
 * so, how to make the theory easy?
 *
 */
using namespace std;

/*
 * simulation and graphical display of coloumbic potential of a Si atom.
 * eqn -> V = (1/4*pi*eps0) * (q/r)
 * consider 1d grid, then 2d and then 3d.
 * 1d grid, treat the atom as at origin -> as all the atoms are the same, for a series of atoms,
 * we don't need to recompute the potential again and again.
 *
 * Other quantities -> distance between nuclie -> 15A.
 * compute the nuclear potential due a nucleus till 5A distance from the nucleus.
 * take steps of 0.1A initially.
 *
 * e=1.609 * e-19 C
 * 1/4*pi*eps0 = 9*e9 N-m^2
 * Use GSL constants for all std const values.
 *
 * Object design -> nothing much. need an grid array of the entire length with 5A displacement
 * need to store the potential at each point.
 *
 *
 */

#include <gsl/gsl_vector.h>
#include <gsl/gsl_const_mksa.h>
#include <fstream>


/**
 * Function declarations
 */
void formPotandSave(gsl_vector* vect, char* filename,bool save);
void formEField(gsl_vector* vect, char* filename, bool save);


/**
 * Global variable and constant declaration
 */
double a = 1e-12;
double aL = 5e-10;
int npoints = aL/a;
const double e = GSL_CONST_MKSA_ELECTRON_CHARGE;
double e_mass = GSL_CONST_MKSA_MASS_ELECTRON;
double eps = GSL_CONST_MKSA_VACUUM_PERMITTIVITY;
const double epsConst = 9e9;


int main(){

	cout<<"Testing"<<endl;
	/*
	 * 1.form the single atom potential over a radial distance of 5A
	 */

	cout<<"The a, the aL and the npoints are "<<a<<" ,"<<aL<<" ,"<<npoints<<endl;

	gsl_vector* potential;
	gsl_vector* efield;
	formPotandSave(potential,"pot.dat",true);
//	formEField(efield,"efield.dat",true);

	/*
	 * 2.apply the above calculated potential over the entire grid length, with 5A gap between the effects of each nucleus.
	 */



	return 0;
}


/**
 * Forms the electrical potential at each point and saves to file
 * if required
 */
void formPotandSave(gsl_vector* vect, char* filename,bool save){
	vect = gsl_vector_alloc(npoints);
	int i = 0;
	double pot = 0;
	double d = a;
	ofstream output;

	if(save){
		output.open(filename);
	}
	for(i=0;i<npoints;i++){
		pot = (epsConst * e) /d;
		gsl_vector_set(vect,i,pot);
		if(save){
			output<<d<<"\t"<<gsl_vector_get(vect,i)<<endl;
		}
		d+=a;
	}

	if(save){
		output.close();
	}

	return;
}

/**
 * forms the electrical field strength at each point and saves to a file if required.
 */
void formEField(gsl_vector* vect, char* filename, bool save){
	vect = gsl_vector_alloc(npoints);
	int i=0;
	double efield = 0;
	double d= a;

	ofstream output;
	if(save){
		output.open(filename);
	}

	for(i=0;i<npoints;i++){
		efield = (epsConst*e*e)/(d*d);
		gsl_vector_set(vect,i,efield);
		if(save){
			output<<d<<"\t"<<gsl_vector_get(vect,i)<<endl;
		}
		d+=a;
	}

	if(save){
		output.close();
	}

}





































